!> author: 左志华
!> date: 2022-07-04
!>
!> Smoothed kernel function <br>
!> 光滑核函数
module smoothed_kernel_function_m

    use math_m, only: rp, r23
    implicit none
    private

    public :: cubic_spline_kernel, scale, selfden
    integer :: scale        !! Scale of the kernel function
                            !! 光滑长度缩尺比例
    real(rp) :: selfden     !! Self-density of the particle
                            !! 自身密度
    
    real(rp), parameter :: Pi = acos(-1.0_rp)   !! Pi <br>
                                                !! 圆周率


contains

    !> Smoothed kernel function: Cubic spline <br>
    !> 四次样条曲线 (monaghan 1985) <br>
    !> 计算光滑函数 \(W_{ij}\) 及其导数 \(\frac{dW}{dx_{ij}}\) 的子例程
    pure subroutine cubic_spline_kernel(r, dx, hsml, w, dwdx)
        real(rp), intent(in) :: r           !! Euclidean distance between two points <br>
                                            !! 欧氏距离
        real(rp), intent(in) :: dx(2)       !! Dimensional distance between two points <br>
                                            !! 对应坐标轴距离
        real(rp), intent(in) :: hsml        !! Smoothing length <br>
                                            !! 光滑长度
        real(rp), intent(out) :: w          !! Kernel value <br>
                                            !! 光滑函数
        real(rp), intent(out) :: dwdx(2)    !! Derivative of kernel value <br>
                                            !! 光滑函数导数

        real(rp) :: factor, q

        q = r/hsml
        factor = 15.0_rp/(7.0_rp*Pi*hsml*hsml)

        if (q <= 1.0_rp) then
            w = factor*(r23 - q*q + q**3/2)
            dwdx = factor*(-2.0_rp + 1.5*q)*dx/(hsml*hsml)     ! 计算机原则，先乘后除
        else if (q > 1.0_rp) then
            w = factor*(2.0_rp - q)**3/6.0_rp
            dwdx = -factor*3.0_rp*(2.0_rp - q)**2*dx/(6.0_rp*hsml*r)
        end if

    end subroutine cubic_spline_kernel

end module smoothed_kernel_function_m
